home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Languguage OS 2
/
Languguage OS II Version 10-94 (Knowledge Media)(1994).ISO
/
language
/
asa
/
asa.c
< prev
next >
Wrap
C/C++ Source or Header
|
1994-10-23
|
93KB
|
3,163 lines
/***********************************************************************
* Lester Ingber (copyright) (c)
* See COPYING License in this directory
* Date 1 Jan 93
***********************************************************************/
#define ASA_ID \
"/* $Id: asa.c,v 4.2 1994/10/23 23:35:16 ingber Exp ingber $ */"
#include "asa.h"
/* static variables in EXPONENT_CHECK () */
static double MIN_EXPONENT, MAX_EXPONENT;
/***********************************************************************
* asa
* This procedure implements the full ASA function optimization.
***********************************************************************/
#if HAVE_ANSI
double
asa (double (*user_cost_function) (
double *, double *, double *, double *, double *,
ALLOC_INT *, int *, int *, int *, USER_DEFINES *),
double (*user_random_generator) (),
double *parameter_initial_final,
double *parameter_minimum,
double *parameter_maximum,
double *tangents,
double *curvature,
ALLOC_INT * number_parameters,
int *parameter_type,
int *valid_state_generated_flag,
int *exit_status,
USER_DEFINES * OPTIONS)
#else
double
asa (user_cost_function,
user_random_generator,
parameter_initial_final,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
OPTIONS)
double (*user_cost_function) ();
double (*user_random_generator) ();
double *parameter_initial_final;
double *parameter_minimum;
double *parameter_maximum;
double *tangents;
double *curvature;
ALLOC_INT *number_parameters;
int *parameter_type;
int *valid_state_generated_flag;
int *exit_status;
USER_DEFINES *OPTIONS;
#endif
{
int index_cost_constraint, /* index initial cost functions averaged
to set initial temperature */
tmp_var_int, /* integer temporary value */
index_cost_repeat; /* test OPTIONS->COST_PRECISION when =
OPTIONS->MAXIMUM_COST_REPEAT */
ALLOC_INT index_v, /* iteration index */
*start_sequence; /* initial OPTIONS->SEQUENTIAL_PARAMETERS
used if >= 0 */
double sampled_cost_sum, /* temporary initial costs */
final_cost, /* best cost to return to user */
tmp_var_db1, tmp_var_db2; /* temporary doubles */
int *curvature_flag;
FILE *ptr_asa_out; /* file ptr to output file */
/* The 3 states that are kept track of during the annealing process */
STATE *current_generated_state, *last_saved_state, *best_generated_state;
#if ASA_PARALLEL
/* parallel generated states */
STATE *gener_block_state;
#endif
/* Holds values of parameter increments */
double *delta_parameter;
/* Holds values of quench_param scales */
double *quench_param;
/* Holds value of quench_cost scales */
double *quench_cost;
/* The array of tangents (absolute value of the numerical derivatives),
and the maximum |tangent| of the array */
double *maximum_tangent;
/* ratio of acceptances to generated points - determines when to
test/reanneal */
double *accepted_to_generated_ratio;
/* temperature parameters */
double temperature_scale, *temperature_scale_parameters;
/* relative scalings of cost and parameters to temperature_scale */
double *temperature_scale_cost;
double *current_user_parameter_temp;
double *initial_user_parameter_temp;
double *current_cost_temperature;
double *initial_cost_temperature;
double log_new_temperature_ratio; /* current *temp = initial *temp *
exp(log_new_temperature_ratio) */
ALLOC_INT *index_exit_v; /* information for asa_exit */
/* counts of generated states and acceptances */
LONG_INT *index_parameter_generations;
LONG_INT *number_generated, *best_number_generated_saved;
LONG_INT *recent_number_generated, *number_accepted;
LONG_INT *recent_number_acceptances, *index_cost_acceptances;
LONG_INT *number_acceptances_saved, *best_number_accepted_saved;
/* Flag indicates that the parameters generated were
invalid according to the cost function validity criteria. */
LONG_INT *number_invalid_generated_states;
LONG_INT repeated_invalid_states;
#if ASA_PARALLEL
LONG_INT index_parallel; /* count of parallel generated states */
LONG_INT parallel_generated; /* saved *recent_number_generated */
LONG_INT parallel_block_max; /* saved OPTIONS->gener_block_max */
#endif
/* used to index repeated and recursive calls to asa */
/* This assumes that multiple calls (>= 1) _or_ recursive
calls are being made to asa */
static int asa_open = FALSE;
static int number_asa_open = 0;
static int recursive_asa_open = 0;
/* static variables used in EXPONENT_CHECK() */
MIN_EXPONENT = 0.9 * log ((double) MIN_DOUBLE);
MAX_EXPONENT = 0.9 * log ((double) MAX_DOUBLE);
/* initializations */
if ((curvature_flag =
(int *) calloc (1, sizeof (int))) == NULL)
exit (9);
if ((maximum_tangent =
(double *) calloc (1, sizeof (double))) == NULL)
exit (9);
if ((accepted_to_generated_ratio =
(double *) calloc (1, sizeof (double))) == NULL)
exit (9);
if ((temperature_scale_cost =
(double *) calloc (1, sizeof (double))) == NULL)
exit (9);
if ((current_cost_temperature =
(double *) calloc (1, sizeof (double))) == NULL)
exit (9);
if ((initial_cost_temperature =
(double *) calloc (1, sizeof (double))) == NULL)
exit (9);
if ((index_exit_v =
(ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL)
exit (9);
if ((start_sequence =
(ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL)
exit (9);
if ((number_generated =
(LONG_INT *) calloc (1, sizeof (LONG_INT))) == NULL)
exit (9);
if ((best_number_generated_saved =
(LONG_INT *) calloc (1, sizeof (LONG_INT))) == NULL)
exit (9);
if ((recent_number_generated =
(LONG_INT *) calloc (1, sizeof (LONG_INT))) == NULL)
exit (9);
if ((number_accepted =
(LONG_INT *) calloc (1, sizeof (LONG_INT))) == NULL)
exit (9);
if ((recent_number_acceptances =
(LONG_INT *) calloc (1, sizeof (LONG_INT))) == NULL)
exit (9);
if ((index_cost_acceptances =
(LONG_INT *) calloc (1, sizeof (LONG_INT))) == NULL)
exit (9);
if ((number_acceptances_saved =
(LONG_INT *) calloc (1, sizeof (LONG_INT))) == NULL)
exit (9);
if ((best_number_accepted_saved =
(LONG_INT *) calloc (1, sizeof (LONG_INT))) == NULL)
exit (9);
if ((number_invalid_generated_states =
(LONG_INT *) calloc (1, sizeof (LONG_INT))) == NULL)
exit (9);
if ((current_generated_state =
(STATE *) calloc (1, sizeof (STATE))) == NULL)
exit (9);
if ((last_saved_state =
(STATE *) calloc (1, sizeof (STATE))) == NULL)
exit (9);
if ((best_generated_state =
(STATE *) calloc (1, sizeof (STATE))) == NULL)
exit (9);
#if ASA_PARALLEL
if ((gener_block_state =
(STATE *) calloc (OPTIONS->gener_block_max, sizeof (STATE))) == NULL)
exit (9);
#endif
if (asa_open == FALSE)
{
asa_open = TRUE;
++number_asa_open;
#if ASA_PRINT
if (number_asa_open == 1)
{
/* open the output file */
#if USER_ASA_OUT
if (!strcmp (OPTIONS->asa_out_file, "STDOUT"))
{
ptr_asa_out = stdout;
}
else
{
ptr_asa_out = fopen (OPTIONS->asa_out_file, "w");
}
#else
if (!strcmp (ASA_OUT, "STDOUT"))
{
ptr_asa_out = stdout;
}
else
{
ptr_asa_out = fopen (ASA_OUT, "w");
}
#endif
}
else
{
#if USER_ASA_OUT
ptr_asa_out = fopen (OPTIONS->asa_out_file, "a");
#else
ptr_asa_out = fopen (ASA_OUT, "a");
#endif
fprintf (ptr_asa_out, "\n\n\t\t number_asa_open = %d\n",
number_asa_open);
}
#endif
}
else
{
++recursive_asa_open;
#if ASA_PRINT
if (recursive_asa_open == 1)
{
/* open the output file */
#if USER_ASA_OUT
ptr_asa_out = fopen (OPTIONS->asa_out_file, "w");
#else
ptr_asa_out = fopen (ASA_OUT, "w");
#endif
}
else
{
#if USER_ASA_OUT
ptr_asa_out = fopen (OPTIONS->asa_out_file, "a");
#else
ptr_asa_out = fopen (ASA_OUT, "a");
#endif
fprintf (ptr_asa_out, "\n\n\t\t recursive_asa_open = %d\n",
recursive_asa_open);
}
#endif
}
#if ASA_PRINT
/* print header information as defined by user */
print_asa_options (ptr_asa_out, OPTIONS);
#if TIME_CALC
/* print starting time */
print_time ("start_asa", ptr_asa_out);
#endif
fflush (ptr_asa_out);
#endif
/* set indices and counts to 0 */
*best_number_generated_saved =
*number_generated =
*recent_number_generated =
*recent_number_acceptances = 0;
*index_cost_acceptances =
*best_number_accepted_saved =
*number_accepted =
*number_acceptances_saved = 0;
index_cost_repeat = 0;
#if ASA_SAMPLE
OPTIONS->n_accepted = 0;
OPTIONS->average_weights = 1.0;
#endif
/* do not calculate curvatures initially */
*curvature_flag = FALSE;
/* allocate storage for all parameters */
if ((current_generated_state->parameter =
(double *) calloc (*number_parameters, sizeof (double))
) == NULL)
exit (9);
if ((last_saved_state->parameter =
(double *) calloc (*number_parameters, sizeof (double))
) == NULL)
exit (9);
if ((best_generated_state->parameter =
(double *) calloc (*number_parameters, sizeof (double))
) == NULL)
exit (9);
#if ASA_PARALLEL
parallel_block_max = OPTIONS->gener_block_max;
parallel_generated = OPTIONS->gener_block;
for (index_parallel = 0; index_parallel < parallel_block_max;
++index_parallel)
{
if ((gener_block_state[index_parallel].parameter =
(double *) calloc (*number_parameters, sizeof (double))
) == NULL)
exit (9);
}
#endif
if ((initial_user_parameter_temp =
(double *) calloc (*number_parameters, sizeof (double))
) == NULL)
exit (9);
if ((index_parameter_generations =
(LONG_INT *) calloc (*number_parameters, sizeof (LONG_INT))
) == NULL)
exit (9);
/* set all delta_parameter[] */
if (OPTIONS->DELTA_PARAMETERS == TRUE)
{
delta_parameter = OPTIONS->user_delta_parameter;
VFOR (index_v)
delta_parameter[index_v] =
OPTIONS->user_delta_parameter[index_v];
}
else
{
if ((delta_parameter =
(double *) calloc (*number_parameters, sizeof (double))
) == NULL)
exit (9);
VFOR (index_v)
delta_parameter[index_v] = OPTIONS->DELTA_X;
}
/* set all quench_param[] */
if (OPTIONS->QUENCH_PARAMETERS == TRUE)
{
quench_param = OPTIONS->user_quench_param_scale;
VFOR (index_v)
quench_param[index_v] =
OPTIONS->user_quench_param_scale[index_v];
}
else
{
if ((quench_param =
(double *) calloc (*number_parameters, sizeof (double))
) == NULL)
exit (9);
VFOR (index_v)
quench_param[index_v] = ONE;
}
/* set quench_cost[] */
if (OPTIONS->QUENCH_COST == TRUE)
{
quench_cost = OPTIONS->user_quench_cost_scale;
*quench_cost = OPTIONS->user_quench_cost_scale[0];
}
else
{
if ((quench_cost =
(double *) calloc (1, sizeof (double))) == NULL)
exit (9);
*quench_cost = ONE;
}
/* set all temperatures */
if ((current_user_parameter_temp =
(double *) calloc (*number_parameters, sizeof (double))
) == NULL)
exit (9);
if (OPTIONS->USER_INITIAL_PARAMETERS_TEMPS == TRUE)
{
VFOR (index_v)
current_user_parameter_temp[index_v] =
initial_user_parameter_temp[index_v] =
OPTIONS->user_parameter_temperature[index_v];
}
else
{
VFOR (index_v)
current_user_parameter_temp[index_v] =
initial_user_parameter_temp[index_v] =
OPTIONS->INITIAL_PARAMETER_TEMPERATURE;
}
if ((temperature_scale_parameters =
(double *) calloc (*number_parameters, sizeof (double))
) == NULL)
exit (9);
if (OPTIONS->USER_INITIAL_COST_TEMP == TRUE)
*initial_cost_temperature = *current_cost_temperature =
OPTIONS->user_cost_temperature[0];
/* set parameters to the initial parameter values */
VFOR (index_v)
last_saved_state->parameter[index_v] =
current_generated_state->parameter[index_v] =
parameter_initial_final[index_v];
/* test OPTIONS->SEQUENTIAL_PARAMETERS */
if (OPTIONS->SEQUENTIAL_PARAMETERS >= *number_parameters)
{
#if ASA_PRINT
fprintf (ptr_asa_out, "OPTIONS->SEQUENTIAL_PARAMETERS > dimension-1\n");
fflush (ptr_asa_out);
#endif
exit (9);
}
/* save initial user value of OPTIONS->SEQUENTIAL_PARAMETERS */
*start_sequence = OPTIONS->SEQUENTIAL_PARAMETERS;
#if ASA_PRINT
#if INT_ALLOC
fprintf (ptr_asa_out,
"*number_parameters = %d\n\n", *number_parameters);
#else
fprintf (ptr_asa_out,
#if INT_LONG
"*number_parameters = %ld\n\n", *number_parameters);
#else
"*number_parameters = %d\n\n", *number_parameters);
#endif
#endif
/* print the min, max, current values, and types of parameters */
fprintf (ptr_asa_out,
"index_v parameter_minimum parameter_maximum\
parameter_value parameter_type \n");
#if ASA_PRINT_INTERMED
VFOR (index_v)
#if INT_ALLOC
fprintf (ptr_asa_out, " %-8d %-18.7g %-17.7g %-15.7g %-7d\n",
#else
#if INT_LONG
fprintf (ptr_asa_out, " %-8ld %-18.7g %-17.7g %-15.7g %-7d\n",
#else
fprintf (ptr_asa_out, " %-8d %-18.7g %-17.7g %-15.7g %-7d\n",
#endif
#endif
index_v,
parameter_minimum[index_v],
parameter_maximum[index_v],
current_generated_state->parameter[index_v],
parameter_type[index_v]);
fprintf (ptr_asa_out, "\n\n");
#endif
/* Print out user-defined OPTIONS */
if (OPTIONS->DELTA_PARAMETERS)
VFOR (index_v)
fprintf (ptr_asa_out,
#if INT_ALLOC
"OPTIONS->user_delta_parameter[%d] = %12.7g\n",
#else
#if INT_LONG
"OPTIONS->user_delta_parameter[%ld] = %12.7g\n",
#else
"OPTIONS->user_delta_parameter[%d] = %12.7g\n",
#endif
#endif
index_v, delta_parameter[index_v]);
if (OPTIONS->QUENCH_PARAMETERS)
VFOR (index_v)
fprintf (ptr_asa_out,
#if INT_ALLOC
"OPTIONS->user_quench_param_scale[%d] = %12.7g\n",
#else
#if INT_LONG
"OPTIONS->user_quench_param_scale[%ld] = %12.7g\n",
#else
"OPTIONS->user_quench_param_scale[%d] = %12.7g\n",
#endif
#endif
index_v, quench_param[index_v]);
if (OPTIONS->QUENCH_COST)
fprintf (ptr_asa_out,
"OPTIONS->user_quench_cost_scale = %12.7g\n",
*quench_cost);
if (OPTIONS->USER_INITIAL_PARAMETERS_TEMPS)
VFOR (index_v)
fprintf (ptr_asa_out,
#if INT_ALLOC
"OPTIONS->user_parameter_temperature[%d] = %12.7g\n",
#else
#if INT_LONG
"OPTIONS->user_parameter_temperature[%ld] = %12.7g\n",
#else
"OPTIONS->user_parameter_temperature[%d] = %12.7g\n",
#endif
#endif
index_v, initial_user_parameter_temp[index_v]);
if (OPTIONS->RATIO_TEMPERATURE_SCALES)
VFOR (index_v)
fprintf (ptr_asa_out,
#if INT_ALLOC
"OPTIONS->user_temperature_ratio[%d] = %12.7g\n",
#else
#if INT_LONG
"OPTIONS->user_temperature_ratio[%ld] = %12.7g\n",
#else
"OPTIONS->user_temperature_ratio[%d] = %12.7g\n",
#endif
#endif
index_v, OPTIONS->user_temperature_ratio[index_v]);
if (OPTIONS->USER_INITIAL_COST_TEMP)
fprintf (ptr_asa_out,
"OPTIONS->user_cost_temperature[0] = %12.7g\n",
*initial_cost_temperature);
fflush (ptr_asa_out);
#endif
#if ASA_PARALLEL
#if ASA_PRINT
fprintf (ptr_asa_out,
#if INT_LONG
"Initial ASA_PARALLEL OPTIONS->\n\t gener_block = %ld\n\
\t gener_block_max = %ld\n \t gener_mov_avr= %d\n\n",
#else
"ASA_PARALLEL OPTIONS->\n\t gener_block = %d\n\
\t gener_block_max = %d\n \t gener_mov_avr= %d\n\n",
#endif
OPTIONS->gener_block, OPTIONS->gener_block_max, OPTIONS->gener_mov_avr);
#endif
#endif
/* set the initial index of parameter generations to 1 */
VFOR (index_v)
index_parameter_generations[index_v] = 1;
if (OPTIONS->USER_INITIAL_COST_TEMP == FALSE)
{
/* calculate the average cost over samplings of the cost function */
sampled_cost_sum = ZERO;
for (index_cost_constraint = 0;
index_cost_constraint < (OPTIONS->NUMBER_COST_SAMPLES);
++index_cost_constraint)
{
*number_invalid_generated_states = 0;
repeated_invalid_states = 0;
OPTIONS->SEQUENTIAL_PARAMETERS = *start_sequence - 1;
do
{
++(*number_invalid_generated_states);
++repeated_invalid_states;
generate_new_state (user_random_generator,
parameter_minimum,
parameter_maximum,
current_user_parameter_temp,
number_parameters,
parameter_type,
current_generated_state,
last_saved_state,
OPTIONS);
*valid_state_generated_flag = TRUE;
current_generated_state->cost =
user_cost_function (current_generated_state->parameter,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
OPTIONS);
if (repeated_invalid_states >
OPTIONS->LIMIT_INVALID_GENERATED_STATES)
{
*exit_status = (int) TOO_MANY_INVALID_STATES;
goto EXIT_ASA;
}
}
while (*valid_state_generated_flag == FALSE);
--(*number_invalid_generated_states);
sampled_cost_sum += fabs (current_generated_state->cost);
}
/* set all cost temperatures to the average cost */
*initial_cost_temperature = *current_cost_temperature =
sampled_cost_sum / (OPTIONS->NUMBER_COST_SAMPLES);
}
/* set all parameters to the initial parameter values */
VFOR (index_v)
best_generated_state->parameter[index_v] =
last_saved_state->parameter[index_v] =
current_generated_state->parameter[index_v] =
parameter_initial_final[index_v];
/* if using user's initial parameters */
if (OPTIONS->USER_INITIAL_PARAMETERS == TRUE)
{
*valid_state_generated_flag = TRUE;
current_generated_state->cost =
user_cost_function (current_generated_state->parameter,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
OPTIONS);
#if ASA_PRINT
if (*valid_state_generated_flag == FALSE)
fprintf (ptr_asa_out,
"user's initial parameters generated \
FALSE *valid_state_generated_flag\n");
#endif
}
else
{
/* let asa generate valid initial parameters */
repeated_invalid_states = 0;
OPTIONS->SEQUENTIAL_PARAMETERS = *start_sequence - 1;
do
{
++(*number_invalid_generated_states);
++repeated_invalid_states;
generate_new_state (user_random_generator,
parameter_minimum,
parameter_maximum,
current_user_parameter_temp,
number_parameters,
parameter_type,
current_generated_state,
last_saved_state,
OPTIONS);
*valid_state_generated_flag = TRUE;
current_generated_state->cost =
user_cost_function (current_generated_state->parameter,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
OPTIONS);
if (repeated_invalid_states >
OPTIONS->LIMIT_INVALID_GENERATED_STATES)
{
*exit_status = (int) TOO_MANY_INVALID_STATES;
goto EXIT_ASA;
}
}
while (*valid_state_generated_flag == FALSE);
--(*number_invalid_generated_states);
} /* OPTIONS->USER_INITIAL_PARAMETERS */
/* set all states to the last one generated */
VFOR (index_v)
{
/* ignore parameters that have too small a range */
if (PARAMETER_RANGE_TOO_SMALL (index_v))
continue;
best_generated_state->parameter[index_v] =
last_saved_state->parameter[index_v] =
current_generated_state->parameter[index_v];
}
/* set all costs to the last one generated */
best_generated_state->cost = last_saved_state->cost =
current_generated_state->cost;
*accepted_to_generated_ratio = ONE;
tmp_var_db1 = -log ((OPTIONS->TEMPERATURE_RATIO_SCALE));
/* compute temperature scales */
tmp_var_db2 = log (OPTIONS->TEMPERATURE_ANNEAL_SCALE);
temperature_scale = tmp_var_db1 * exp (-tmp_var_db2
/ (double) *number_parameters);
if (OPTIONS->RATIO_TEMPERATURE_SCALES)
{
VFOR (index_v)
temperature_scale_parameters[index_v] =
tmp_var_db1 * exp (-(tmp_var_db2 * quench_param[index_v])
/ (double) *number_parameters)
* OPTIONS->user_temperature_ratio[index_v];
}
else
{
VFOR (index_v)
temperature_scale_parameters[index_v] =
tmp_var_db1 * exp (-(tmp_var_db2 * quench_param[index_v])
/ (double) *number_parameters);
}
*temperature_scale_cost =
tmp_var_db1 * exp (-(tmp_var_db2 * quench_cost[0])
/ (double) *number_parameters)
* OPTIONS->COST_PARAMETER_SCALE;
/* do not calculate curvatures initially */
*curvature_flag = FALSE;
#if ASA_PRINT
fprintf (ptr_asa_out,
"temperature_scale = %12.7g\n",
temperature_scale);
if (OPTIONS->RATIO_TEMPERATURE_SCALES)
{
#if ASA_PRINT_INTERMED
VFOR (index_v)
{
fprintf (ptr_asa_out,
#if INT_ALLOC
"temperature_scale_parameters[%d] = %12.7g\n",
#else
#if INT_LONG
"temperature_scale_parameters[%ld] = %12.7g\n",
#else
"temperature_scale_parameters[%d] = %12.7g\n",
#endif
#endif
index_v, temperature_scale_parameters[index_v]);
}
#endif
}
else
{
fprintf (ptr_asa_out,
"temperature_scale_parameters[0] = %12.7g\n",
temperature_scale_parameters[0]);
}
fprintf (ptr_asa_out,
"*temperature_scale_cost = %12.7g\n",
*temperature_scale_cost);
fprintf (ptr_asa_out, "\n\n");
#if ASA_PRINT_INTERMED
print_state (parameter_minimum,
parameter_maximum,
tangents,
curvature,
current_cost_temperature,
current_user_parameter_temp,
accepted_to_generated_ratio,
number_parameters,
curvature_flag,
number_accepted,
index_cost_acceptances,
number_generated,
number_invalid_generated_states,
last_saved_state,
best_generated_state,
ptr_asa_out,
OPTIONS);
#endif
fprintf (ptr_asa_out, "\n");
fflush (ptr_asa_out);
#endif
#if ASA_SAMPLE
#if ASA_PRINT
fprintf (ptr_asa_out,
":SAMPLE: n_accept cost cost_temp bias_accept \
aver_weight\n");
fprintf (ptr_asa_out,
":SAMPLE: index param[] temp[] bias_gener[] \
range[]\n");
#endif
#endif
/* reset the current cost and the number of generations performed */
*number_invalid_generated_states = 0;
*best_number_generated_saved =
*number_generated = *recent_number_generated = 0;
VFOR (index_v)
{
/* ignore parameters that have too small a range */
if (PARAMETER_RANGE_TOO_SMALL (index_v))
continue;
index_parameter_generations[index_v] = 1;
}
OPTIONS->SEQUENTIAL_PARAMETERS = *start_sequence - 1;
/* MAIN ANNEALING LOOP */
while (
((*number_accepted <= OPTIONS->LIMIT_ACCEPTANCES)
|| (OPTIONS->LIMIT_ACCEPTANCES == 0))
&&
((*number_generated <= OPTIONS->LIMIT_GENERATED)
|| (OPTIONS->LIMIT_GENERATED == 0))
)
{
tmp_var_db1 = -log ((OPTIONS->TEMPERATURE_RATIO_SCALE));
/* compute temperature scales */
tmp_var_db2 = log (OPTIONS->TEMPERATURE_ANNEAL_SCALE);
temperature_scale = tmp_var_db1 * exp (-tmp_var_db2
/ (double) *number_parameters);
if (OPTIONS->RATIO_TEMPERATURE_SCALES)
{
VFOR (index_v)
temperature_scale_parameters[index_v] =
tmp_var_db1 * exp (-(tmp_var_db2 * quench_param[index_v])
/ (double) *number_parameters)
* OPTIONS->user_temperature_ratio[index_v];
}
else
{
VFOR (index_v)
temperature_scale_parameters[index_v] =
tmp_var_db1 * exp (-(tmp_var_db2 * quench_param[index_v])
/ (double) *number_parameters);
}
*temperature_scale_cost =
tmp_var_db1 * exp (-(tmp_var_db2 * quench_cost[0])
/ (double) *number_parameters)
* OPTIONS->COST_PARAMETER_SCALE;
/* CALCULATE NEW TEMPERATURES */
/* calculate new parameter temperatures */
VFOR (index_v)
{
/* skip parameters with too small a range */
if (PARAMETER_RANGE_TOO_SMALL (index_v))
continue;
log_new_temperature_ratio =
-temperature_scale_parameters[index_v] *
pow ((double) index_parameter_generations[index_v],
quench_param[index_v] / (double) *number_parameters);
/* check (and correct) for too large an exponent */
log_new_temperature_ratio =
EXPONENT_CHECK (log_new_temperature_ratio);
current_user_parameter_temp[index_v] =
initial_user_parameter_temp[index_v]
* exp (log_new_temperature_ratio);
#if NO_PARAM_TEMP_TEST
if (current_user_parameter_temp[index_v]
< (double) EPS_DOUBLE)
current_user_parameter_temp[index_v] =
(double) EPS_DOUBLE;
#else
/* check for too small a parameter temperature */
if (current_user_parameter_temp[index_v]
< (double) EPS_DOUBLE)
{
*exit_status = (int) P_TEMP_TOO_SMALL;
*index_exit_v = index_v;
goto EXIT_ASA;
}
#endif
}
/* calculate new cost temperature */
log_new_temperature_ratio =
-*temperature_scale_cost
* pow ((double) *index_cost_acceptances,
(*quench_cost) / (double) *number_parameters);
log_new_temperature_ratio =
EXPONENT_CHECK (log_new_temperature_ratio);
*current_cost_temperature = *initial_cost_temperature
* exp (log_new_temperature_ratio);
#if NO_COST_TEMP_TEST
if (*current_cost_temperature < (double) EPS_DOUBLE)
*current_cost_temperature = (double) EPS_DOUBLE;
#else
/* check for too small a cost temperature */
if (*current_cost_temperature < (double) EPS_DOUBLE)
{
*exit_status = (int) C_TEMP_TOO_SMALL;
goto EXIT_ASA;
}
#endif
/* GENERATE NEW PARAMETERS */
/* generate a new valid set of parameters */
#if ASA_PARALLEL
/* *** ENTER CODE TO SPAWN OFF PARALLEL GENERATED STATES *** */
if (OPTIONS->gener_block_max > parallel_block_max)
{
for (index_parallel = 0; index_parallel < parallel_block_max;
++index_parallel)
{
free (gener_block_state[index_parallel].parameter);
}
free (gener_block_state);
if ((gener_block_state =
(STATE *) calloc (OPTIONS->gener_block_max,
sizeof (STATE))) == NULL)
exit (9);
parallel_block_max = OPTIONS->gener_block_max;
for (index_parallel = 0; index_parallel < parallel_block_max;
++index_parallel)
{
if ((gener_block_state[index_parallel].parameter =
(double *) calloc (*number_parameters, sizeof (double))
) == NULL)
exit (9);
}
}
for (index_parallel = 0; index_parallel < OPTIONS->gener_block;
++index_parallel)
{
#endif /* ASA_PARALLEL */
repeated_invalid_states = 0;
do
{
++(*number_invalid_generated_states);
++repeated_invalid_states;
generate_new_state (user_random_generator,
parameter_minimum,
parameter_maximum,
current_user_parameter_temp,
number_parameters,
parameter_type,
current_generated_state,
last_saved_state,
OPTIONS);
*valid_state_generated_flag = TRUE;
current_generated_state->cost =
user_cost_function (current_generated_state->parameter,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
OPTIONS);
if (repeated_invalid_states >
OPTIONS->LIMIT_INVALID_GENERATED_STATES)
{
*exit_status = (int) TOO_MANY_INVALID_STATES;
goto EXIT_ASA;
}
}
while (*valid_state_generated_flag == FALSE);
--(*number_invalid_generated_states);
#if ASA_PARALLEL
gener_block_state[index_parallel].cost
= current_generated_state->cost;
VFOR (index_v)
{
/* ignore parameters with too small a range */
if (PARAMETER_RANGE_TOO_SMALL (index_v))
continue;
gener_block_state[index_parallel].parameter[index_v] =
current_generated_state->parameter[index_v];
}
}
/* *** EXIT CODE SPAWNING OFF PARALLEL GENERATED STATES *** */
#endif /* ASA_PARALLEL */
/* ACCEPT/REJECT NEW PARAMETERS */
#if ASA_PARALLEL
for (index_parallel = 0; index_parallel < OPTIONS->gener_block;
++index_parallel)
{
current_generated_state->cost
= gener_block_state[index_parallel].cost;
VFOR (index_v)
{
/* ignore parameters with too small a range */
if (PARAMETER_RANGE_TOO_SMALL (index_v))
continue;
current_generated_state->parameter[index_v] =
gener_block_state[index_parallel].parameter[index_v];
}
#endif /* ASA_PARALLEL */
/* decide to accept/reject the new state */
accept_new_state (user_random_generator,
parameter_minimum,
parameter_maximum,
current_cost_temperature,
#if ASA_SAMPLE
current_user_parameter_temp,
#endif
number_parameters,
recent_number_acceptances,
number_accepted,
index_cost_acceptances,
number_acceptances_saved,
recent_number_generated,
number_generated,
index_parameter_generations,
current_generated_state,
last_saved_state,
#if ASA_SAMPLE
ptr_asa_out,
#endif
OPTIONS);
/* calculate the ratio of acceptances to generated states */
*accepted_to_generated_ratio =
(double) (*recent_number_acceptances + 1) /
(double) (*recent_number_generated + 1);
/* CHECK FOR NEW MINIMUM */
if (current_generated_state->cost < best_generated_state->cost)
{
/* NEW MINIMUM FOUND */
/* reset the recent acceptances and generated counts */
#if ASA_PARALLEL
parallel_generated = *recent_number_generated;
#endif
*recent_number_acceptances = *recent_number_generated = 0;
*best_number_generated_saved = *number_generated;
*best_number_accepted_saved = *number_accepted;
index_cost_repeat = 0;
/* copy the current state into the best_generated state */
best_generated_state->cost = current_generated_state->cost;
VFOR (index_v)
{
/* ignore parameters that have too small a range */
if (PARAMETER_RANGE_TOO_SMALL (index_v))
continue;
best_generated_state->parameter[index_v] =
current_generated_state->parameter[index_v];
}
/* printout the new minimum state and value */
#if ASA_PRINT
#if INT_LONG
fprintf (ptr_asa_out,
"best...->cost=%-12.7g \
*number_accepted=%ld *number_generated=%ld\n",
best_generated_state->cost,
*number_accepted, *number_generated);
#else
fprintf (ptr_asa_out,
"best...->cost=%-12.7g \
*number_accepted=%d *number_generated=%d\n",
best_generated_state->cost,
*number_accepted, *number_generated);
#endif /* INT_LOG */
#if ASA_PARALLEL
/* print OPTIONS->gener_block just used */
fprintf (ptr_asa_out,
#if INT_LONG
"OPTIONS->gener_block = %ld\n", OPTIONS->gener_block);
#else
"OPTIONS->gener_block = %d\n", OPTIONS->gener_block);
#endif /* INT_LONG */
#endif /* ASA_PARALLEL */
#if ASA_PRINT_MORE
print_state (parameter_minimum,
parameter_maximum,
tangents,
curvature,
current_cost_temperature,
current_user_parameter_temp,
accepted_to_generated_ratio,
number_parameters,
curvature_flag,
number_accepted,
index_cost_acceptances,
number_generated,
number_invalid_generated_states,
last_saved_state,
best_generated_state,
ptr_asa_out,
OPTIONS);
#endif /* ASA_PRINT_MORE */
fflush (ptr_asa_out);
#endif /* ASA_PRINT */
#if ASA_PARALLEL
/* leave index_parallel loop after new minimum */
break;
#endif /* ASA_PARALLEL */
}
#if ASA_PARALLEL
}
#endif /* ASA_PARALLEL */
#if ASA_PARALLEL
if (OPTIONS->gener_mov_avr > 0)
{
OPTIONS->gener_block = (LONG_INT)
((((double) OPTIONS->gener_mov_avr - ONE)
* (double) (OPTIONS->gener_block) + (double) parallel_generated)
/ (double) (OPTIONS->gener_mov_avr));
OPTIONS->gener_block =
MIN (OPTIONS->gener_block, parallel_block_max);
}
#endif /* ASA_PARALLEL */
/* PERIODIC TESTING/REANNEALING/PRINTING SECTION */
tmp_var_int = (int) (*number_accepted %
((LONG_INT) OPTIONS->TESTING_FREQUENCY_MODULUS));
if ((tmp_var_int == 0 && *number_acceptances_saved
== *number_accepted)
|| *accepted_to_generated_ratio
< (OPTIONS->ACCEPTED_TO_GENERATED_RATIO))
{
if (*accepted_to_generated_ratio
< (OPTIONS->ACCEPTED_TO_GENERATED_RATIO))
*recent_number_acceptances = *recent_number_generated = 0;
/* if best.cost repeats OPTIONS->MAXIMUM_COST_REPEAT then exit */
if (fabs (last_saved_state->cost
- best_generated_state->cost) <
OPTIONS->COST_PRECISION)
{
++index_cost_repeat;
if (index_cost_repeat == (OPTIONS->MAXIMUM_COST_REPEAT))
{
*exit_status = (int) COST_REPEATING;
goto EXIT_ASA;
}
}
else
{
index_cost_repeat = 0;
}
/* set to FALSE to skip reannealing */
if (OPTIONS->ACTIVATE_REANNEAL == TRUE)
{
/* calculate tangents, not curvatures, to reanneal */
*curvature_flag = FALSE;
cost_derivatives (user_cost_function,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
maximum_tangent,
delta_parameter,
number_parameters,
parameter_type,
exit_status,
curvature_flag,
valid_state_generated_flag,
number_invalid_generated_states,
current_generated_state,
best_generated_state,
ptr_asa_out,
OPTIONS);
reanneal (parameter_minimum,
parameter_maximum,
tangents,
maximum_tangent,
current_cost_temperature,
initial_cost_temperature,
temperature_scale_cost,
current_user_parameter_temp,
initial_user_parameter_temp,
temperature_scale_parameters,
quench_param,
quench_cost,
number_parameters,
parameter_type,
index_cost_acceptances,
index_parameter_generations,
last_saved_state,
best_generated_state,
OPTIONS);
}
#if ASA_PRINT_INTERMED
#if ASA_PRINT
print_state (parameter_minimum,
parameter_maximum,
tangents,
curvature,
current_cost_temperature,
current_user_parameter_temp,
accepted_to_generated_ratio,
number_parameters,
curvature_flag,
number_accepted,
index_cost_acceptances,
number_generated,
number_invalid_generated_states,
last_saved_state,
best_generated_state,
ptr_asa_out,
OPTIONS);
fprintf (ptr_asa_out, "\n");
fflush (ptr_asa_out);
#endif
#endif
}
}
/* FINISHED ANNEALING and MINIMIZATION */
*exit_status = (int) NORMAL_EXIT;
EXIT_ASA:
asa_exit (user_cost_function,
&final_cost,
parameter_initial_final,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
maximum_tangent,
delta_parameter,
current_cost_temperature,
initial_user_parameter_temp,
current_user_parameter_temp,
accepted_to_generated_ratio,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
index_exit_v,
start_sequence,
number_accepted,
best_number_accepted_saved,
index_cost_acceptances,
number_generated,
number_invalid_generated_states,
index_parameter_generations,
best_number_generated_saved,
current_generated_state,
last_saved_state,
best_generated_state,
ptr_asa_out,
OPTIONS);
free (curvature_flag);
free (maximum_tangent);
free (accepted_to_generated_ratio);
free (temperature_scale_cost);
free (current_cost_temperature);
free (initial_cost_temperature);
free (number_generated);
free (best_number_generated_saved);
free (recent_number_generated);
free (number_accepted);
free (recent_number_acceptances);
free (index_cost_acceptances);
free (number_acceptances_saved);
free (best_number_accepted_saved);
free (number_invalid_generated_states);
free (current_generated_state->parameter);
free (last_saved_state->parameter);
free (best_generated_state->parameter);
free (current_generated_state);
free (last_saved_state);
free (best_generated_state);
#if ASA_PARALLEL
for (index_parallel = 0; index_parallel < parallel_block_max;
++index_parallel)
{
free (gener_block_state[index_parallel].parameter);
}
free (gener_block_state);
#endif
free (initial_user_parameter_temp);
free (index_exit_v);
free (start_sequence);
free (index_parameter_generations);
if (OPTIONS->DELTA_PARAMETERS == FALSE)
free (delta_parameter);
if (OPTIONS->QUENCH_PARAMETERS == FALSE)
free (quench_param);
if (OPTIONS->QUENCH_COST == FALSE)
free (quench_cost);
free (current_user_parameter_temp);
free (temperature_scale_parameters);
if (recursive_asa_open == 0)
asa_open = FALSE;
return (final_cost);
}
/***********************************************************************
* asa_exit
* This procedures copies the best parameters and cost into
* final_cost and parameter_initial_final
***********************************************************************/
#if HAVE_ANSI
void
asa_exit (double (*user_cost_function) (
double *, double *, double *, double *, double *,
ALLOC_INT *, int *, int *, int *, USER_DEFINES *),
double *final_cost,
double *parameter_initial_final,
double *parameter_minimum,
double *parameter_maximum,
double *tangents,
double *curvature,
double *maximum_tangent,
double *delta_parameter,
double *current_cost_temperature,
double *initial_user_parameter_temp,
double *current_user_parameter_temp,
double *accepted_to_generated_ratio,
ALLOC_INT * number_parameters,
int *parameter_type,
int *valid_state_generated_flag,
int *exit_status,
ALLOC_INT * index_exit_v,
ALLOC_INT * start_sequence,
LONG_INT * number_accepted,
LONG_INT * best_number_accepted_saved,
LONG_INT * index_cost_acceptances,
LONG_INT * number_generated,
LONG_INT * number_invalid_generated_states,
LONG_INT * index_parameter_generations,
LONG_INT * best_number_generated_saved,
STATE * current_generated_state,
STATE * last_saved_state,
STATE * best_generated_state,
FILE * ptr_asa_out,
USER_DEFINES * OPTIONS)
#else
void
asa_exit (user_cost_function,
final_cost,
parameter_initial_final,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
maximum_tangent,
delta_parameter,
current_cost_temperature,
initial_user_parameter_temp,
current_user_parameter_temp,
accepted_to_generated_ratio,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
index_exit_v,
start_sequence,
number_accepted,
best_number_accepted_saved,
index_cost_acceptances,
number_generated,
number_invalid_generated_states,
index_parameter_generations,
best_number_generated_saved,
current_generated_state,
last_saved_state,
best_generated_state,
ptr_asa_out,
OPTIONS)
double (*user_cost_function) ();
double *final_cost;
double *parameter_initial_final;
double *parameter_minimum;
double *parameter_maximum;
double *tangents;
double *curvature;
double *maximum_tangent;
double *delta_parameter;
double *current_cost_temperature;
double *initial_user_parameter_temp;
double *current_user_parameter_temp;
double *accepted_to_generated_ratio;
ALLOC_INT *number_parameters;
int *parameter_type;
int *valid_state_generated_flag;
int *exit_status;
ALLOC_INT *index_exit_v;
ALLOC_INT *start_sequence;
LONG_INT *number_accepted;
LONG_INT *best_number_accepted_saved;
LONG_INT *index_cost_acceptances;
LONG_INT *number_generated;
LONG_INT *number_invalid_generated_states;
LONG_INT *index_parameter_generations;
LONG_INT *best_number_generated_saved;
STATE *current_generated_state;
STATE *last_saved_state;
STATE *best_generated_state;
FILE *ptr_asa_out;
USER_DEFINES *OPTIONS;
#endif
{
ALLOC_INT index_v; /* iteration index */
int *curvature_flag;
if ((curvature_flag =
(int *) calloc (1, sizeof (int))) == NULL)
exit (9);
if (*exit_status != TOO_MANY_INVALID_STATES)
{
/* calculate curvatures and tangents at best point */
*curvature_flag = TRUE;
cost_derivatives (user_cost_function,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
maximum_tangent,
delta_parameter,
number_parameters,
parameter_type,
exit_status,
curvature_flag,
valid_state_generated_flag,
number_invalid_generated_states,
current_generated_state,
best_generated_state,
ptr_asa_out,
OPTIONS);
}
/* return final function minimum and associated parameters */
*final_cost = best_generated_state->cost;
VFOR (index_v)
{
parameter_initial_final[index_v] =
best_generated_state->parameter[index_v];
}
#if ASA_PRINT
print_state (parameter_minimum,
parameter_maximum,
tangents,
curvature,
current_cost_temperature,
current_user_parameter_temp,
accepted_to_generated_ratio,
number_parameters,
curvature_flag,
number_accepted,
index_cost_acceptances,
number_generated,
number_invalid_generated_states,
last_saved_state,
best_generated_state,
ptr_asa_out,
OPTIONS);
switch (*exit_status)
{
case NORMAL_EXIT:
fprintf (ptr_asa_out,
"\n\n NORMAL_EXIT exit_status = %d\n",
*exit_status);
break;
case P_TEMP_TOO_SMALL:
fprintf (ptr_asa_out,
"\n\n P_TEMP_TOO_SMALL exit_status = %d\n",
*exit_status);
fprintf (ptr_asa_out,
#if INT_ALLOC
"current_user_parameter_temp[%d] too small = %12.7g\n",
#else
#if INT_LONG
"current_user_parameter_temp[%ld] too small = %12.7g\n",
#else
"current_user_parameter_temp[%d] too small = %12.7g\n",
#endif
#endif
*index_exit_v, current_user_parameter_temp[*index_exit_v]);
break;
case C_TEMP_TOO_SMALL:
fprintf (ptr_asa_out,
"\n\n C_TEMP_TOO_SMALL exit_status = %d\n",
*exit_status);
fprintf (ptr_asa_out,
"*current_cost_temperature too small = %12.7g\n",
*current_cost_temperature);
break;
case COST_REPEATING:
fprintf (ptr_asa_out,
"\n\n COST_REPEATING exit_status = %d\n",
*exit_status);
break;
case TOO_MANY_INVALID_STATES:
fprintf (ptr_asa_out,
"\n\n TOO_MANY_INVALID_STATES exit_status = %d\n",
*exit_status);
break;
default:
fprintf (ptr_asa_out, "\n\n ERR: no exit code available = %d\n",
*exit_status);
}
fprintf (ptr_asa_out,
"final_cost = best_generated_state->cost = %-12.7g\n",
*final_cost);
#if INT_LONG
fprintf (ptr_asa_out,
"*number_accepted at best_generated_state->cost = %ld\n",
*best_number_accepted_saved);
fprintf (ptr_asa_out,
"*number_generated at best_generated_state->cost = %ld\n",
*best_number_generated_saved);
#else
fprintf (ptr_asa_out,
"*number_accepted at best_generated_state->cost = %d\n",
*best_number_accepted_saved);
fprintf (ptr_asa_out,
"*number_generated at best_generated_state->cost = %d\n",
*best_number_generated_saved);
#endif
#endif
#if ASA_TEMPLATE
#if OPTIONAL_DATA
if (OPTIONS->asa_data[0] > (double) MIN_DOUBLE)
OPTIONS->asa_data[1] = (double) (*best_number_generated_saved);
#endif
#endif
/* reset OPTIONS->SEQUENTIAL_PARAMETERS */
OPTIONS->SEQUENTIAL_PARAMETERS = *start_sequence;
/* return unused space */
free (curvature_flag);
#if ASA_PRINT
#if TIME_CALC
/* print ending time */
print_time ("asa_end", ptr_asa_out);
#endif
fprintf (ptr_asa_out, "\n\n\n");
fflush (ptr_asa_out);
fclose (ptr_asa_out);
#endif
}
/***********************************************************************
* generate_new_state
* Generates a valid new state from the old state
***********************************************************************/
#if HAVE_ANSI
void
generate_new_state (double (*user_random_generator) (),
double *parameter_minimum,
double *parameter_maximum,
double *current_user_parameter_temp,
ALLOC_INT * number_parameters,
int *parameter_type,
STATE * current_generated_state,
STATE * last_saved_state,
USER_DEFINES * OPTIONS)
#else
void
generate_new_state (user_random_generator,
parameter_minimum,
parameter_maximum,
current_user_parameter_temp,
number_parameters,
parameter_type,
current_generated_state,
last_saved_state,
OPTIONS)
double (*user_random_generator) ();
double *parameter_minimum;
double *parameter_maximum;
double *current_user_parameter_temp;
ALLOC_INT *number_parameters;
int *parameter_type;
STATE *current_generated_state;
STATE *last_saved_state;
USER_DEFINES *OPTIONS;
#endif
{
ALLOC_INT index_v;
double x;
double parameter_v, min_parameter_v, max_parameter_v, temperature_v,
parameter_range_v;
/* generate a new value for each parameter */
VFOR (index_v)
{
if (OPTIONS->SEQUENTIAL_PARAMETERS >= -1)
{
++OPTIONS->SEQUENTIAL_PARAMETERS;
if (OPTIONS->SEQUENTIAL_PARAMETERS == *number_parameters)
OPTIONS->SEQUENTIAL_PARAMETERS = 0;
index_v = OPTIONS->SEQUENTIAL_PARAMETERS;
}
min_parameter_v = parameter_minimum[index_v];
max_parameter_v = parameter_maximum[index_v];
parameter_range_v = max_parameter_v - min_parameter_v;
/* ignore parameters that have too small a range */
if (fabs (parameter_range_v) < (double) EPS_DOUBLE)
continue;
temperature_v = current_user_parameter_temp[index_v];
parameter_v = last_saved_state->parameter[index_v];
/* Handle discrete integer parameters. */
if (INTEGER_PARAMETER (index_v))
{
min_parameter_v -= HALF;
max_parameter_v += HALF;
parameter_range_v = max_parameter_v - min_parameter_v;
}
/* generate a new state x within the parameter bounds */
for (;;)
{
x = parameter_v
+ generate_asa_state (user_random_generator, &temperature_v)
* parameter_range_v;
/* exit the loop if within its valid parameter range */
if (x <= max_parameter_v - (double) EPS_DOUBLE
&& x >= min_parameter_v + (double) EPS_DOUBLE)
break;
}
/* Handle discrete integer parameters.
You might have to check rounding on your machine. */
if (INTEGER_PARAMETER (index_v))
{
if (x < min_parameter_v + HALF)
x = min_parameter_v + HALF + (double) EPS_DOUBLE;
if (x > max_parameter_v - HALF)
x = max_parameter_v - HALF + (double) EPS_DOUBLE;
if (x + HALF > ZERO)
{
x = (int) (x + HALF);
}
else
{
x = (int) (x - HALF);
}
if (x > parameter_maximum[index_v])
x = parameter_maximum[index_v];
if (x < parameter_minimum[index_v])
x = parameter_minimum[index_v];
}
/* save the newly generated value */
current_generated_state->parameter[index_v] = x;
if (OPTIONS->SEQUENTIAL_PARAMETERS >= 0)
break;
}
}
/***********************************************************************
* generate_asa_state
* This function generates a single value according to the
* ASA generating function and the passed temperature
***********************************************************************/
#if HAVE_ANSI
double
generate_asa_state (double (*user_random_generator) (),
double *temp)
#else
double
generate_asa_state (user_random_generator,
temp)
double (*user_random_generator) ();
double *temp;
#endif
{
double x, y, z;
x = (*user_random_generator) ();
y = x < HALF ? -ONE : ONE;
z = y * *temp * (pow ((ONE + ONE / *temp), fabs (TWO * x - ONE)) - ONE);
return (z);
}
/***********************************************************************
* accept_new_state
* This procedure accepts or rejects a newly generated state,
* depending on whether the difference between new and old
* cost functions passes a statistical test. If accepted,
* the current state is updated.
***********************************************************************/
#if HAVE_ANSI
void
accept_new_state (double (*user_random_generator) (),
double *parameter_minimum,
double *parameter_maximum,
double *current_cost_temperature,
#if ASA_SAMPLE
double *current_user_parameter_temp,
#endif
ALLOC_INT * number_parameters,
LONG_INT * recent_number_acceptances,
LONG_INT * number_accepted,
LONG_INT * index_cost_acceptances,
LONG_INT * number_acceptances_saved,
LONG_INT * recent_number_generated,
LONG_INT * number_generated,
LONG_INT * index_parameter_generations,
STATE * current_generated_state,
STATE * last_saved_state,
#if ASA_SAMPLE
FILE * ptr_asa_out,
#endif
USER_DEFINES * OPTIONS)
#else
void
accept_new_state (user_random_generator,
parameter_minimum,
parameter_maximum,
current_cost_temperature,
#if ASA_SAMPLE
current_user_parameter_temp,
#endif
number_parameters,
recent_number_acceptances,
number_accepted,
index_cost_acceptances,
number_acceptances_saved,
recent_number_generated,
number_generated,
index_parameter_generations,
current_generated_state,
last_saved_state,
#if ASA_SAMPLE
ptr_asa_out,
#endif
OPTIONS)
double (*user_random_generator) ();
double *parameter_minimum;
double *parameter_maximum;
double *current_cost_temperature;
#if ASA_SAMPLE
double *current_user_parameter_temp;
#endif
ALLOC_INT *number_parameters;
LONG_INT *recent_number_acceptances;
LONG_INT *number_accepted;
LONG_INT *index_cost_acceptances;
LONG_INT *number_acceptances_saved;
LONG_INT *recent_number_generated;
LONG_INT *number_generated;
LONG_INT *index_parameter_generations;
STATE *current_generated_state;
STATE *last_saved_state;
#if ASA_SAMPLE
FILE *ptr_asa_out;
#endif
USER_DEFINES *OPTIONS;
#endif
{
double delta_cost, boltz_test, unif_test;
ALLOC_INT index_v;
double curr_cost_temp;
#if ASA_SAMPLE
LONG_INT active_params;
double weight_param_ind, weight_aver, range;
#endif
/* update accepted and generated count */
++*number_acceptances_saved;
++*recent_number_generated;
++*number_generated;
/* increment the parameter index generation for each parameter */
if (OPTIONS->SEQUENTIAL_PARAMETERS >= 0)
{
/* ignore parameters with too small a range */
if (!PARAMETER_RANGE_TOO_SMALL (OPTIONS->SEQUENTIAL_PARAMETERS))
++index_parameter_generations[OPTIONS->SEQUENTIAL_PARAMETERS];
}
else
{
VFOR (index_v)
{
if (!PARAMETER_RANGE_TOO_SMALL (index_v))
++index_parameter_generations[index_v];
}
}
/* effective cost function for testing acceptance criteria,
calculate the cost difference and divide by the temperature */
#if USER_COST_SCHEDULE
curr_cost_temp =
(OPTIONS->cost_schedule (*current_cost_temperature, OPTIONS)
+ EPS_DOUBLE);
#else
curr_cost_temp = *current_cost_temperature;
#endif
delta_cost = (current_generated_state->cost - last_saved_state->cost)
/ (curr_cost_temp + EPS_DOUBLE);
boltz_test = MIN (ONE, (exp (EXPONENT_CHECK (-delta_cost))));
unif_test = (*user_random_generator) ();
#if ASA_SAMPLE
active_params = 0;
weight_aver = ZERO;
VFOR (index_v)
{
/* ignore parameters with too small a range */
if (PARAMETER_RANGE_TOO_SMALL (index_v))
continue;
++active_params;
range = parameter_maximum[index_v] - parameter_minimum[index_v];
weight_param_ind = TWO
* (fabs ((last_saved_state->parameter[index_v]
- current_generated_state->parameter[index_v]) / (TWO * range))
+ current_user_parameter_temp[index_v])
* log (ONE + ONE / current_user_parameter_temp[index_v]);
weight_aver += weight_param_ind;
OPTIONS->bias_generated[index_v] = ONE / weight_param_ind;
}
weight_aver /= (double) active_params;
OPTIONS->average_weights = weight_aver;
OPTIONS->bias_acceptance = boltz_test;
if (boltz_test >= unif_test)
{
++OPTIONS->n_accepted;
}
#if ASA_PRINT
if (OPTIONS->limit_weights < OPTIONS->average_weights)
{
fprintf (ptr_asa_out, ":SAMPLE#\n");
if (boltz_test >= unif_test)
{
#if INT_LONG
fprintf (ptr_asa_out, ":SAMPLE+ %10ld %12.7g %12.7g %12.7g %12.7g\n",
#else
fprintf (ptr_asa_out, ":SAMPLE+ %10d %12.7g %12.7g %12.7g\n",
#endif
OPTIONS->n_accepted, current_generated_state->cost,
curr_cost_temp, OPTIONS->bias_acceptance,
OPTIONS->average_weights);
VFOR (index_v)
{
/* ignore parameters with too small a range */
if (PARAMETER_RANGE_TOO_SMALL (index_v))
continue;
range = parameter_maximum[index_v] - parameter_minimum[index_v];
#if INT_ALLOC
fprintf (ptr_asa_out, ":SAMPLE %11d %12.7g %12.7g %12.7g %12.7g\n",
#else
#if INT_LONG
fprintf (ptr_asa_out, ":SAMPLE %11ld %12.7g %12.7g %12.7g %12.7g\n",
#else
fprintf (ptr_asa_out, ":SAMPLE %11d %12.7g %12.7g %12.7g %12.7g\n",
#endif
#endif
index_v, current_generated_state->parameter[index_v],
current_user_parameter_temp[index_v],
OPTIONS->bias_generated[index_v], range);
}
}
else
{
#if INT_LONG
fprintf (ptr_asa_out, ":SAMPLE %11ld %12.7g %12.7g %12.7g %12.7g\n",
#else
fprintf (ptr_asa_out, ":SAMPLE %11d %12.7g %12.7g %12.7g\n",
#endif
OPTIONS->n_accepted, last_saved_state->cost,
curr_cost_temp, OPTIONS->bias_acceptance,
OPTIONS->average_weights);
VFOR (index_v)
{
/* ignore parameters with too small a range */
if (PARAMETER_RANGE_TOO_SMALL (index_v))
continue;
range = parameter_maximum[index_v] - parameter_minimum[index_v];
#if INT_ALLOC
fprintf (ptr_asa_out, ":SAMPLE %11d %12.7g %12.7g %12.7g %12.7g\n",
#else
#if INT_LONG
fprintf (ptr_asa_out, ":SAMPLE %11ld %12.7g %12.7g %12.7g %12.7g\n",
#else
fprintf (ptr_asa_out, ":SAMPLE %11d %12.7g %12.7g %12.7g %12.7g\n",
#endif
#endif
index_v, last_saved_state->parameter[index_v],
current_user_parameter_temp[index_v],
OPTIONS->bias_generated[index_v], range);
}
}
}
#endif
#endif /* ASA_SAMPLE */
/* accept/reject the new state */
if (boltz_test >= unif_test)
{
/* copy current state to the last saved state */
last_saved_state->cost = current_generated_state->cost;
VFOR (index_v)
{
/* ignore parameters with too small a range */
if (PARAMETER_RANGE_TOO_SMALL (index_v))
continue;
last_saved_state->parameter[index_v] =
current_generated_state->parameter[index_v];
}
/* update acceptance counts */
++*recent_number_acceptances;
++*number_accepted;
++*index_cost_acceptances;
*number_acceptances_saved = *number_accepted;
}
}
/***********************************************************************
* reanneal
* Readjust temperatures of generating and acceptance functions
***********************************************************************/
#if HAVE_ANSI
void
reanneal (double *parameter_minimum,
double *parameter_maximum,
double *tangents,
double *maximum_tangent,
double *current_cost_temperature,
double *initial_cost_temperature,
double *temperature_scale_cost,
double *current_user_parameter_temp,
double *initial_user_parameter_temp,
double *temperature_scale_parameters,
double *quench_param,
double *quench_cost,
ALLOC_INT * number_parameters,
int *parameter_type,
LONG_INT * index_cost_acceptances,
LONG_INT * index_parameter_generations,
STATE * last_saved_state,
STATE * best_generated_state,
USER_DEFINES * OPTIONS)
#else
void
reanneal (parameter_minimum,
parameter_maximum,
tangents,
maximum_tangent,
current_cost_temperature,
initial_cost_temperature,
temperature_scale_cost,
current_user_parameter_temp,
initial_user_parameter_temp,
temperature_scale_parameters,
quench_param,
quench_cost,
number_parameters,
parameter_type,
index_cost_acceptances,
index_parameter_generations,
last_saved_state,
best_generated_state,
OPTIONS)
double *parameter_minimum;
double *parameter_maximum;
double *tangents;
double *maximum_tangent;
double *current_cost_temperature;
double *initial_cost_temperature;
double *temperature_scale_cost;
double *current_user_parameter_temp;
double *initial_user_parameter_temp;
double *temperature_scale_parameters;
double *quench_param;
double *quench_cost;
ALLOC_INT *number_parameters;
int *parameter_type;
LONG_INT *index_cost_acceptances;
LONG_INT *index_parameter_generations;
STATE *last_saved_state;
STATE *best_generated_state;
USER_DEFINES *OPTIONS;
#endif
{
ALLOC_INT index_v;
double tmp_var_db3;
double new_temperature; /* the temperature */
double log_new_temperature_ratio;
double log_init_cur_temp_ratio;
double temperature_rescale_power;
double tmp_dbl;
VFOR (index_v)
{
if (NO_REANNEAL (index_v))
continue;
/* use the temp double to prevent overflow */
tmp_dbl = (double) index_parameter_generations[index_v];
/* skip parameters with too small range or integer parameters */
if (OPTIONS->INCLUDE_INTEGER_PARAMETERS == TRUE)
{
if (PARAMETER_RANGE_TOO_SMALL (index_v))
continue;
}
else
{
if (PARAMETER_RANGE_TOO_SMALL (index_v) ||
INTEGER_PARAMETER (index_v))
continue;
}
/* ignore parameters with too small tangents */
if (fabs (tangents[index_v]) < (double) EPS_DOUBLE)
continue;
/* reset the index of parameter generations appropriately */
#if USER_REANNEAL_FUNCTION
new_temperature =
fabs (OPTIONS->reanneal_function (current_user_parameter_temp[index_v],
tangents[index_v],
*maximum_tangent,
OPTIONS));
#else
new_temperature =
fabs (REANNEAL_FUNCTION (current_user_parameter_temp[index_v],
tangents[index_v],
*maximum_tangent));
#endif
if (new_temperature < initial_user_parameter_temp[index_v])
{
log_init_cur_temp_ratio =
fabs (log (((double) EPS_DOUBLE
+ initial_user_parameter_temp[index_v])
/ ((double) EPS_DOUBLE + new_temperature)));
tmp_dbl = (double) EPS_DOUBLE
+ pow (log_init_cur_temp_ratio
/ temperature_scale_parameters[index_v],
(double) *number_parameters
/ quench_param[index_v]);
}
else
{
tmp_dbl = ONE;
}
/* Reset index_parameter_generations if index reset too large,
and also reset the initial_user_parameter_temp, to achieve
the same new temperature. */
while (tmp_dbl
> ((double) OPTIONS->MAXIMUM_REANNEAL_INDEX))
{
log_new_temperature_ratio =
-temperature_scale_parameters[index_v] *
pow (tmp_dbl,
quench_param[index_v] / (double) *number_parameters);
log_new_temperature_ratio =
EXPONENT_CHECK (log_new_temperature_ratio);
new_temperature =
initial_user_parameter_temp[index_v]
* exp (log_new_temperature_ratio);
tmp_dbl /=
OPTIONS->REANNEAL_RESCALE;
temperature_rescale_power = ONE
/ pow (OPTIONS->REANNEAL_RESCALE,
quench_param[index_v]
/ (double) *number_parameters);
initial_user_parameter_temp[index_v] =
new_temperature * pow (initial_user_parameter_temp
[index_v] / new_temperature,
temperature_rescale_power);
}
/* restore from temporary double */
index_parameter_generations[index_v] = (LONG_INT) tmp_dbl;
}
/* reanneal : Reset the index of cost acceptances to take
into account finer detail in cost terrain. */
/* set the starting cost_temperature to last cost found so far */
if (*initial_cost_temperature > fabs (last_saved_state->cost))
*initial_cost_temperature = fabs (last_saved_state->cost);
tmp_dbl = (double) *index_cost_acceptances;
if (*current_cost_temperature > fabs (best_generated_state->cost))
{
tmp_var_db3 =
fabs (log (((double) EPS_DOUBLE + *initial_cost_temperature) /
((double) EPS_DOUBLE + fabs (best_generated_state->cost))));
tmp_dbl = (double) EPS_DOUBLE
+ pow (tmp_var_db3
/ *temperature_scale_cost,
(double) *number_parameters / (*quench_cost));
}
else
{
log_init_cur_temp_ratio =
fabs (log (((double) EPS_DOUBLE + *initial_cost_temperature) /
((double) EPS_DOUBLE + *current_cost_temperature)));
tmp_dbl = (double) EPS_DOUBLE
+ pow (log_init_cur_temp_ratio
/ *temperature_scale_cost,
(double) *number_parameters / (*quench_cost));
}
/* reset index_cost_temperature if index reset too large */
while (tmp_dbl > ((double) OPTIONS->MAXIMUM_REANNEAL_INDEX))
{
log_new_temperature_ratio = -*temperature_scale_cost *
pow ((double) tmp_dbl,
(*quench_cost) / (double) *number_parameters);
log_new_temperature_ratio =
EXPONENT_CHECK (log_new_temperature_ratio);
new_temperature = *initial_cost_temperature
* exp (log_new_temperature_ratio);
tmp_dbl /= OPTIONS->REANNEAL_RESCALE;
temperature_rescale_power = ONE
/ pow (OPTIONS->REANNEAL_RESCALE,
(*quench_cost)
/ (double) *number_parameters);
*initial_cost_temperature =
new_temperature * pow (*initial_cost_temperature
/ new_temperature,
temperature_rescale_power);
}
*index_cost_acceptances = (LONG_INT) tmp_dbl;
}
/***********************************************************************
* cost_derivatives
* This procedure calculates the derivatives of the cost function
* with respect to its parameters. The first derivatives are
* used as a sensitivity measure for reannealing. The second
* derivatives are calculated only if *curvature_flag=TRUE;
* these are a measure of the covariance of the fit when a
* minimum is found.
***********************************************************************/
/* Calculate the numerical derivatives of the best
generated state found so far */
/* In this implementation of ASA, no checks are made for
*valid_state_generated_flag=FALSE for differential neighbors
to the current best state. */
/* Assuming no information is given about the metric of the parameter
space, use simple Cartesian space to calculate curvatures. */
#if HAVE_ANSI
void
cost_derivatives (double (*user_cost_function) (
double *, double *, double *, double *, double *,
ALLOC_INT *, int *, int *, int *, USER_DEFINES *),
double *parameter_minimum,
double *parameter_maximum,
double *tangents,
double *curvature,
double *maximum_tangent,
double *delta_parameter,
ALLOC_INT * number_parameters,
int *parameter_type,
int *exit_status,
int *curvature_flag,
int *valid_state_generated_flag,
LONG_INT * number_invalid_generated_states,
STATE * current_generated_state,
STATE * best_generated_state,
FILE * ptr_asa_out,
USER_DEFINES * OPTIONS)
#else
void
cost_derivatives (user_cost_function,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
maximum_tangent,
delta_parameter,
number_parameters,
parameter_type,
exit_status,
curvature_flag,
valid_state_generated_flag,
number_invalid_generated_states,
current_generated_state,
best_generated_state,
ptr_asa_out,
OPTIONS)
double (*user_cost_function) ();
double *parameter_minimum;
double *parameter_maximum;
double *tangents;
double *curvature;
double *maximum_tangent;
double *delta_parameter;
ALLOC_INT *number_parameters;
int *parameter_type;
int *exit_status;
int *curvature_flag;
int *valid_state_generated_flag;
LONG_INT *number_invalid_generated_states;
STATE *current_generated_state;
STATE *best_generated_state;
FILE *ptr_asa_out;
USER_DEFINES *OPTIONS;
#endif
{
ALLOC_INT index_v, index_vv, index_v_vv, index_vv_v;
LONG_INT saved_num_invalid_gen_states, tmp_saved;
double parameter_v, parameter_vv, parameter_v_offset, parameter_vv_offset;
double recent_best_cost;
double new_cost_state_1, new_cost_state_2, new_cost_state_3;
double delta_parameter_v, delta_parameter_vv;
if (OPTIONS->CURVATURE_0 == TRUE)
*curvature_flag = FALSE;
if (OPTIONS->CURVATURE_0 == -1)
*curvature_flag = TRUE;
/* save the best cost */
recent_best_cost = best_generated_state->cost;
/* copy the best state into the current state */
VFOR (index_v)
{
/* ignore parameters with too small ranges */
if (PARAMETER_RANGE_TOO_SMALL (index_v))
continue;
current_generated_state->parameter[index_v] =
best_generated_state->parameter[index_v];
}
saved_num_invalid_gen_states = (*number_invalid_generated_states);
/* set parameters (& possibly constraints) to best state */
*valid_state_generated_flag = TRUE;
current_generated_state->cost =
user_cost_function (current_generated_state->parameter,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
OPTIONS);
if (*valid_state_generated_flag == FALSE)
++(*number_invalid_generated_states);
if (OPTIONS->USER_TANGENTS == TRUE)
{
*valid_state_generated_flag = FALSE;
current_generated_state->cost =
user_cost_function (current_generated_state->parameter,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
OPTIONS);
if (*valid_state_generated_flag == FALSE)
++(*number_invalid_generated_states);
}
else
{
/* calculate tangents */
VFOR (index_v)
{
if (NO_REANNEAL (index_v))
{
tangents[index_v] = ZERO;
continue;
}
/* skip parameters with too small range or integer parameters */
if (OPTIONS->INCLUDE_INTEGER_PARAMETERS == TRUE)
{
if (PARAMETER_RANGE_TOO_SMALL (index_v))
{
tangents[index_v] = ZERO;
continue;
}
}
else
{
if (PARAMETER_RANGE_TOO_SMALL (index_v) ||
INTEGER_PARAMETER (index_v))
{
tangents[index_v] = ZERO;
continue;
}
}
/* save the v_th parameter and delta_parameter */
parameter_v = best_generated_state->parameter[index_v];
delta_parameter_v = delta_parameter[index_v];
parameter_v_offset = (ONE + delta_parameter_v) * parameter_v;
if (parameter_v_offset > parameter_maximum[index_v] ||
parameter_v_offset < parameter_minimum[index_v])
{
delta_parameter_v = -delta_parameter_v;
parameter_v_offset = (ONE + delta_parameter_v) * parameter_v;
}
/* generate the first sample point */
current_generated_state->parameter[index_v] = parameter_v_offset;
*valid_state_generated_flag = TRUE;
current_generated_state->cost =
user_cost_function (current_generated_state->parameter,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
OPTIONS);
if (*valid_state_generated_flag == FALSE)
++(*number_invalid_generated_states);
new_cost_state_1 = current_generated_state->cost;
/* restore the parameter state */
current_generated_state->parameter[index_v] = parameter_v;
/* calculate the numerical derivative */
tangents[index_v] = (new_cost_state_1 - recent_best_cost)
/ (delta_parameter_v * parameter_v + (double) EPS_DOUBLE);
}
}
/* find the maximum |tangent| from all tangents */
*maximum_tangent = 0;
VFOR (index_v)
{
if (NO_REANNEAL (index_v))
continue;
/* ignore too small ranges and integer parameters types */
if (OPTIONS->INCLUDE_INTEGER_PARAMETERS == TRUE)
{
if (PARAMETER_RANGE_TOO_SMALL (index_v))
continue;
}
else
{
if (PARAMETER_RANGE_TOO_SMALL (index_v) ||
INTEGER_PARAMETER (index_v))
continue;
}
/* find the maximum |tangent| (from all tangents) */
if (fabs (tangents[index_v]) > *maximum_tangent)
{
*maximum_tangent = fabs (tangents[index_v]);
}
}
if (*curvature_flag == TRUE || *curvature_flag == -1)
{
/* calculate diagonal curvatures */
VFOR (index_v)
{
if (NO_REANNEAL (index_v))
{
index_v_vv = ROW_COL_INDEX (index_v, index_v);
curvature[index_v_vv] = ZERO;
continue;
}
/* skip parameters with too small range or integer parameters */
if (OPTIONS->INCLUDE_INTEGER_PARAMETERS == TRUE)
{
if (PARAMETER_RANGE_TOO_SMALL (index_v))
{
index_v_vv = ROW_COL_INDEX (index_v, index_v);
curvature[index_v_vv] = ZERO;
continue;
}
else
{
if (PARAMETER_RANGE_TOO_SMALL (index_v) ||
INTEGER_PARAMETER (index_v))
{
index_v_vv = ROW_COL_INDEX (index_v, index_v);
curvature[index_v_vv] = ZERO;
continue;
}
}
}
/* save the v_th parameter and delta_parameter */
parameter_v = best_generated_state->parameter[index_v];
delta_parameter_v = delta_parameter[index_v];
if (parameter_v + delta_parameter_v * fabs (parameter_v)
> parameter_maximum[index_v])
{
/* generate the first sample point */
current_generated_state->parameter[index_v] =
parameter_v - TWO * delta_parameter_v * fabs (parameter_v);
*valid_state_generated_flag = TRUE;
current_generated_state->cost =
user_cost_function (current_generated_state->parameter,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
OPTIONS);
if (*valid_state_generated_flag == FALSE)
++(*number_invalid_generated_states);
new_cost_state_1 = current_generated_state->cost;
/* generate the second sample point */
current_generated_state->parameter[index_v] =
parameter_v - delta_parameter_v * fabs (parameter_v);
*valid_state_generated_flag = TRUE;
current_generated_state->cost =
user_cost_function (current_generated_state->parameter,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
OPTIONS);
if (*valid_state_generated_flag == FALSE)
++(*number_invalid_generated_states);
new_cost_state_2 = current_generated_state->cost;
/* restore the parameter state */
current_generated_state->parameter[index_v] =
parameter_v;
/* index_v_vv: row index_v, column index_v */
index_v_vv = ROW_COL_INDEX (index_v, index_v);
/* calculate and store the curvature */
curvature[index_v_vv] =
(recent_best_cost - TWO * new_cost_state_2
+ new_cost_state_1) / (delta_parameter_v * delta_parameter_v
* parameter_v * parameter_v + (double) EPS_DOUBLE);
}
else if (parameter_v - delta_parameter_v * fabs (parameter_v)
< parameter_minimum[index_v])
{
/* generate the first sample point */
current_generated_state->parameter[index_v] =
parameter_v + TWO * delta_parameter_v * fabs (parameter_v);
*valid_state_generated_flag = TRUE;
current_generated_state->cost =
user_cost_function (current_generated_state->parameter,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
OPTIONS);
if (*valid_state_generated_flag == FALSE)
++(*number_invalid_generated_states);
new_cost_state_1 = current_generated_state->cost;
/* generate the second sample point */
current_generated_state->parameter[index_v] =
parameter_v + delta_parameter_v * fabs (parameter_v);
*valid_state_generated_flag = TRUE;
current_generated_state->cost =
user_cost_function (current_generated_state->parameter,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
OPTIONS);
if (*valid_state_generated_flag == FALSE)
++(*number_invalid_generated_states);
new_cost_state_2 = current_generated_state->cost;
/* restore the parameter state */
current_generated_state->parameter[index_v] =
parameter_v;
/* index_v_vv: row index_v, column index_v */
index_v_vv = ROW_COL_INDEX (index_v, index_v);
/* calculate and store the curvature */
curvature[index_v_vv] =
(recent_best_cost - TWO * new_cost_state_2
+ new_cost_state_1) / (delta_parameter_v * delta_parameter_v
* parameter_v * parameter_v + (double) EPS_DOUBLE);
}
else
{
/* generate the first sample point */
parameter_v_offset = (ONE + delta_parameter_v) * parameter_v;
current_generated_state->parameter[index_v] = parameter_v_offset;
*valid_state_generated_flag = TRUE;
current_generated_state->cost =
user_cost_function (current_generated_state->parameter,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
OPTIONS);
if (*valid_state_generated_flag == FALSE)
++(*number_invalid_generated_states);
new_cost_state_1 = current_generated_state->cost;
/* generate the second sample point */
current_generated_state->parameter[index_v] =
(ONE - delta_parameter_v) * parameter_v;
*valid_state_generated_flag = TRUE;
current_generated_state->cost =
user_cost_function (current_generated_state->parameter,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
OPTIONS);
if (*valid_state_generated_flag == FALSE)
++(*number_invalid_generated_states);
new_cost_state_2 = current_generated_state->cost;
/* restore the parameter state */
current_generated_state->parameter[index_v] =
parameter_v;
/* index_v_vv: row index_v, column index_v */
index_v_vv = ROW_COL_INDEX (index_v, index_v);
/* calculate and store the curvature */
curvature[index_v_vv] =
(new_cost_state_2 - TWO * recent_best_cost
+ new_cost_state_1) / (delta_parameter_v * delta_parameter_v
* parameter_v * parameter_v + (double) EPS_DOUBLE);
}
}
/* calculate off-diagonal curvatures */
VFOR (index_v)
{
/* save the v_th parameter and delta_x */
parameter_v = current_generated_state->parameter[index_v];
delta_parameter_v = delta_parameter[index_v];
VFOR (index_vv)
{
/* index_v_vv: row index_v, column index_vv */
index_v_vv = ROW_COL_INDEX (index_v, index_vv);
index_vv_v = ROW_COL_INDEX (index_vv, index_v);
if (NO_REANNEAL (index_vv) || NO_REANNEAL (index_v))
{
curvature[index_vv_v] = curvature[index_v_vv] = ZERO;
continue;
}
/* calculate only the upper diagonal */
if (index_v <= index_vv)
continue;
/* skip parms with too small range or integer parameters */
if (OPTIONS->INCLUDE_INTEGER_PARAMETERS == TRUE)
{
if (PARAMETER_RANGE_TOO_SMALL (index_v) ||
PARAMETER_RANGE_TOO_SMALL (index_vv))
{
curvature[index_vv_v] = curvature[index_v_vv] = ZERO;
continue;
}
}
else
{
if (INTEGER_PARAMETER (index_v) ||
INTEGER_PARAMETER (index_vv) ||
PARAMETER_RANGE_TOO_SMALL (index_v) ||
PARAMETER_RANGE_TOO_SMALL (index_vv))
{
curvature[index_vv_v] = curvature[index_v_vv] = ZERO;
continue;
}
}
/* save the vv_th parameter and delta_parameter */
parameter_vv =
current_generated_state->parameter[index_vv];
delta_parameter_vv = delta_parameter[index_vv];
/* generate first sample point */
parameter_v_offset = current_generated_state->parameter[index_v] =
(ONE + delta_parameter_v) * parameter_v;
parameter_vv_offset = current_generated_state->parameter[index_vv] =
(ONE + delta_parameter_vv) * parameter_vv;
if (parameter_v_offset > parameter_maximum[index_v] ||
parameter_v_offset < parameter_minimum[index_v])
{
delta_parameter_v = -delta_parameter_v;
current_generated_state->parameter[index_v] =
(ONE + delta_parameter_v) * parameter_v;
}
if (parameter_vv_offset > parameter_maximum[index_vv] ||
parameter_vv_offset < parameter_minimum[index_vv])
{
delta_parameter_vv = -delta_parameter_vv;
current_generated_state->parameter[index_vv] =
(ONE + delta_parameter_vv) * parameter_vv;
}
*valid_state_generated_flag = TRUE;
current_generated_state->cost =
user_cost_function (current_generated_state->parameter,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
OPTIONS);
if (*valid_state_generated_flag == FALSE)
++(*number_invalid_generated_states);
new_cost_state_1 = current_generated_state->cost;
/* restore the v_th parameter */
current_generated_state->parameter[index_v] =
parameter_v;
/* generate second sample point */
*valid_state_generated_flag = TRUE;
current_generated_state->cost =
user_cost_function (current_generated_state->parameter,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
OPTIONS);
if (*valid_state_generated_flag == FALSE)
++(*number_invalid_generated_states);
new_cost_state_2 = current_generated_state->cost;
/* restore the vv_th parameter */
current_generated_state->parameter[index_vv] =
parameter_vv;
/* generate third sample point */
current_generated_state->parameter[index_v] =
(ONE + delta_parameter_v) * parameter_v;
*valid_state_generated_flag = TRUE;
current_generated_state->cost =
user_cost_function (current_generated_state->parameter,
parameter_minimum,
parameter_maximum,
tangents,
curvature,
number_parameters,
parameter_type,
valid_state_generated_flag,
exit_status,
OPTIONS);
if (*valid_state_generated_flag == FALSE)
++(*number_invalid_generated_states);
new_cost_state_3 = current_generated_state->cost;
/* restore the v_th parameter */
current_generated_state->parameter[index_v] =
parameter_v;
/* calculate and store the curvature */
curvature[index_vv_v] = curvature[index_v_vv] =
(new_cost_state_1 - new_cost_state_2
- new_cost_state_3 + recent_best_cost)
/ (delta_parameter_v * delta_parameter_vv
* parameter_v * parameter_vv
+ (double) EPS_DOUBLE);
}
}
}
/* restore the best cost function value */
current_generated_state->cost = recent_best_cost;
#if ASA_PRINT
tmp_saved = *number_invalid_generated_states - saved_num_invalid_gen_states;
if (tmp_saved > 0)
#if INT_LONG
fprintf (ptr_asa_out,
"Generated %ld invalid states when calculating the derivatives\n",
tmp_saved);
#else
fprintf (ptr_asa_out,
"Generated %d invalid states when calculating the derivatives\n",
tmp_saved);
#endif
#endif /* ASA_PRINT */
*number_invalid_generated_states = saved_num_invalid_gen_states;
}
#if ASA_PRINT
/***********************************************************************
* print_state
* Prints a description of the current state of the system
***********************************************************************/
#if HAVE_ANSI
void
print_state (double *parameter_minimum,
double *parameter_maximum,
double *tangents,
double *curvature,
double *current_cost_temperature,
double *current_user_parameter_temp,
double *accepted_to_generated_ratio,
ALLOC_INT * number_parameters,
int *curvature_flag,
LONG_INT * number_accepted,
LONG_INT * index_cost_acceptances,
LONG_INT * number_generated,
LONG_INT * number_invalid_generated_states,
STATE * last_saved_state,
STATE * best_generated_state,
FILE * ptr_asa_out,
USER_DEFINES * OPTIONS)
#else
void
print_state (parameter_minimum,
parameter_maximum,
tangents,
curvature,
current_cost_temperature,
current_user_parameter_temp,
accepted_to_generated_ratio,
number_parameters,
curvature_flag,
number_accepted,
index_cost_acceptances,
number_generated,
number_invalid_generated_states,
last_saved_state,
best_generated_state,
ptr_asa_out,
OPTIONS)
double *parameter_minimum;
double *parameter_maximum;
double *tangents;
double *curvature;
double *current_cost_temperature;
double *current_user_parameter_temp;
double *accepted_to_generated_ratio;
ALLOC_INT *number_parameters;
int *curvature_flag;
LONG_INT *number_accepted;
LONG_INT *index_cost_acceptances;
LONG_INT *number_generated;
LONG_INT *number_invalid_generated_states;
STATE *last_saved_state;
STATE *best_generated_state;
FILE *ptr_asa_out;
USER_DEFINES *OPTIONS;
#endif /* HAVE_ANSI */
{
ALLOC_INT index_v;
ALLOC_INT index_vv, index_v_vv;
fprintf (ptr_asa_out, "\n");
#if TIME_CALC
print_time ("", ptr_asa_out);
#endif
if (OPTIONS->CURVATURE_0 == TRUE)
*curvature_flag = FALSE;
if (OPTIONS->CURVATURE_0 == -1)
*curvature_flag = TRUE;
#if INT_LONG
fprintf (ptr_asa_out,
"*index_cost_acceptances = %ld, *current_cost_temperature = %12.7g\n",
*index_cost_acceptances, *current_cost_temperature);
fprintf (ptr_asa_out,
"*accepted_to_generated_ratio = %12.7g,\
*number_invalid... = %ld\n",
*accepted_to_generated_ratio,
(*number_invalid_generated_states));
fprintf (ptr_asa_out,
"*number_generated = %ld, *number_accepted = %ld\n",
*number_generated, *number_accepted);
#else
fprintf (ptr_asa_out,
"*index_cost_acceptances = %d, *current_cost_temperature = %12.7g\n",
*index_cost_acceptances, *current_cost_temperature);
fprintf (ptr_asa_out,
"*accepted_to_generated_ratio = %12.7g,\
*number_invalid... = %d\n",
*accepted_to_generated_ratio,
*number_invalid_generated_states);
fprintf (ptr_asa_out,
"*number_generated = %d, *number_accepted = %d\n",
*number_generated, *number_accepted);
#endif
fprintf (ptr_asa_out,
"best...->cost = %12.7g,\
last...->cost = %12.7g\n",
best_generated_state->cost, last_saved_state->cost);
/* Note that tangents will not be calculated until reanneal
is called, and therefore their listing in the printout only
is relevant then */
VFOR (index_v)
{
/* ignore too small ranges */
if (PARAMETER_RANGE_TOO_SMALL (index_v))
continue;
fprintf (ptr_asa_out,
#if INT_ALLOC
"best_generated_state->parameter[%d] = %12.7g\n\
current_user_parameter_temp[%d] = %12.7g\n\
tangents[%d]: %12.7g\n",
#else
#if INT_LONG
"best_generated_state->parameter[%ld] = %12.7g\n\
current_user_parameter_temp[%ld] = %12.7g\n\
tangents[%ld]: %12.7g\n",
#else
"best_generated_state->parameter[%d] = %12.7g\n\
current_user_parameter_temp[%d] = %12.7g\n\
tangents[%d]: %12.7g\n",
#endif
#endif
index_v, best_generated_state->parameter[index_v],
index_v, current_user_parameter_temp[index_v],
index_v, tangents[index_v]);
}
if (*curvature_flag == TRUE)
{
/* print curvatures */
VFOR (index_v)
{
/* ignore too small ranges */
if (PARAMETER_RANGE_TOO_SMALL (index_v))
continue;
VFOR (index_vv)
{
/* only print upper diagonal of matrix */
if (index_v < index_vv)
continue;
/* ignore too small ranges (index_vv) */
if (PARAMETER_RANGE_TOO_SMALL (index_vv))
continue;
/* index_v_vv: row index_v, column index_vv */
index_v_vv = ROW_COL_INDEX (index_v, index_vv);
#if INT_ALLOC
fprintf (ptr_asa_out, "curvature[%d][%d] = %12.7g\n",
#else
#if INT_LONG
fprintf (ptr_asa_out, "curvature[%ld][%ld] = %12.7g\n",
#else
fprintf (ptr_asa_out, "curvature[%d][%d] = %12.7g\n",
#endif
#endif
index_v, index_vv, curvature[index_v_vv]);
}
}
}
fprintf (ptr_asa_out, "\n");
fflush (ptr_asa_out);
}
/***********************************************************************
* print_asa-options
* Prints user's selected options
***********************************************************************/
#if HAVE_ANSI
void
print_asa_options (FILE * ptr_asa_out, USER_DEFINES * OPTIONS)
#else
void
print_asa_options (ptr_asa_out, OPTIONS)
FILE *ptr_asa_out;
USER_DEFINES *OPTIONS;
#endif /* HAVE_ANSI */
{
fprintf (ptr_asa_out,
"\t\tADAPTIVE SIMULATED ANNEALING\n\n");
fprintf (ptr_asa_out, "%s\n\n", ASA_ID);
fprintf (ptr_asa_out, "OPTIONS_FILE = %d\n",
(int) OPTIONS_FILE);
fprintf (ptr_asa_out, "ASA_LIB = %d\n",
(int) ASA_LIB);
fprintf (ptr_asa_out, "HAVE_ANSI = %d\n",
(int) HAVE_ANSI);
fprintf (ptr_asa_out, "IO_PROTOTYPES = %d\n",
(int) IO_PROTOTYPES);
fprintf (ptr_asa_out, "TIME_CALC = %d\n",
(int) TIME_CALC);
fprintf (ptr_asa_out, "TIME_STD = %d\n",
(int) TIME_STD);
fprintf (ptr_asa_out, "INT_LONG = %d\n",
(int) INT_LONG);
fprintf (ptr_asa_out, "INT_ALLOC = %d\n",
(int) INT_ALLOC);
fprintf (ptr_asa_out, "SMALL_FLOAT = %g\n",
(double) SMALL_FLOAT);
fprintf (ptr_asa_out, "MIN_DOUBLE = %g\n",
(double) MIN_DOUBLE);
fprintf (ptr_asa_out, "MAX_DOUBLE = %g\n",
(double) MAX_DOUBLE);
fprintf (ptr_asa_out, "EPS_DOUBLE = %g\n",
(double) EPS_DOUBLE);
fprintf (ptr_asa_out, "NO_PARAM_TEMP_TEST = %d\n",
(int) NO_PARAM_TEMP_TEST);
fprintf (ptr_asa_out, "NO_COST_TEMP_TEST = %d\n",
(int) NO_COST_TEMP_TEST);
fprintf (ptr_asa_out, "SELF_OPTIMIZE = %d\n",
(int) SELF_OPTIMIZE);
fprintf (ptr_asa_out, "ASA_TEST = %d\n",
(int) ASA_TEST);
fprintf (ptr_asa_out, "ASA_TEMPLATE = %d\n",
(int) ASA_TEMPLATE);
fprintf (ptr_asa_out, "OPTIONAL_DATA = %d\n",
(int) OPTIONAL_DATA);
fprintf (ptr_asa_out, "USER_COST_SCHEDULE = %d\n",
(int) USER_COST_SCHEDULE);
fprintf (ptr_asa_out, "USER_REANNEAL_FUNCTION = %d\n",
(int) USER_REANNEAL_FUNCTION);
fprintf (ptr_asa_out, "ASA_SAMPLE = %d\n",
(int) ASA_SAMPLE);
fprintf (ptr_asa_out, "ASA_PARALLEL = %d\n\n",
(int) ASA_PARALLEL);
fprintf (ptr_asa_out, "ASA_PRINT = %d\n",
(int) ASA_PRINT);
fprintf (ptr_asa_out, "ASA_OUT = %s\n",
#if USER_ASA_OUT
OPTIONS->asa_out_file);
#else
ASA_OUT);
#endif
fprintf (ptr_asa_out, "USER_ASA_OUT = %d\n",
(int) USER_ASA_OUT);
fprintf (ptr_asa_out, "ASA_PRINT_INTERMED = %d\n",
(int) ASA_PRINT_INTERMED);
fprintf (ptr_asa_out, "ASA_PRINT_MORE = %d\n\n",
(int) ASA_PRINT_MORE);
#if INT_LONG
fprintf (ptr_asa_out, "OPTIONS->LIMIT_ACCEPTANCES = %ld\n",
(LONG_INT) OPTIONS->LIMIT_ACCEPTANCES);
fprintf (ptr_asa_out, "OPTIONS->LIMIT_GENERATED = %ld\n",
(LONG_INT) OPTIONS->LIMIT_GENERATED);
#else
fprintf (ptr_asa_out, "OPTIONS->LIMIT_ACCEPTANCES = %d\n",
(LONG_INT) OPTIONS->LIMIT_ACCEPTANCES);
fprintf (ptr_asa_out, "OPTIONS->LIMIT_GENERATED = %d\n",
(LONG_INT) OPTIONS->LIMIT_GENERATED);
#endif
fprintf (ptr_asa_out, "OPTIONS->LIMIT_INVALID_GENERATED_STATES = %d\n",
OPTIONS->LIMIT_INVALID_GENERATED_STATES);
fprintf (ptr_asa_out, "OPTIONS->ACCEPTED_TO_GENERATED_RATIO = %g\n\n",
OPTIONS->ACCEPTED_TO_GENERATED_RATIO);
fprintf (ptr_asa_out, "OPTIONS->COST_PRECISION = %g\n",
OPTIONS->COST_PRECISION);
fprintf (ptr_asa_out, "OPTIONS->MAXIMUM_COST_REPEAT = %d\n",
OPTIONS->MAXIMUM_COST_REPEAT);
fprintf (ptr_asa_out, "OPTIONS->NUMBER_COST_SAMPLES = %d\n",
OPTIONS->NUMBER_COST_SAMPLES);
fprintf (ptr_asa_out, "OPTIONS->TEMPERATURE_RATIO_SCALE = %g\n",
OPTIONS->TEMPERATURE_RATIO_SCALE);
fprintf (ptr_asa_out, "OPTIONS->COST_PARAMETER_SCALE = %g\n",
OPTIONS->COST_PARAMETER_SCALE);
fprintf (ptr_asa_out, "OPTIONS->TEMPERATURE_ANNEAL_SCALE = %g\n",
OPTIONS->TEMPERATURE_ANNEAL_SCALE);
fprintf (ptr_asa_out, "OPTIONS->USER_INITIAL_COST_TEMP = %d\n\n",
OPTIONS->USER_INITIAL_COST_TEMP);
fprintf (ptr_asa_out, "OPTIONS->INCLUDE_INTEGER_PARAMETERS = %d\n",
OPTIONS->INCLUDE_INTEGER_PARAMETERS);
fprintf (ptr_asa_out, "OPTIONS->USER_INITIAL_PARAMETERS = %d\n",
OPTIONS->USER_INITIAL_PARAMETERS);
#if INT_ALLOC
fprintf (ptr_asa_out, "OPTIONS->SEQUENTIAL_PARAMETERS = %d\n",
(int) OPTIONS->SEQUENTIAL_PARAMETERS);
#else
#if INT_LONG
fprintf (ptr_asa_out, "OPTIONS->SEQUENTIAL_PARAMETERS = %ld\n",
(LONG_INT) OPTIONS->SEQUENTIAL_PARAMETERS);
#else
fprintf (ptr_asa_out, "OPTIONS->SEQUENTIAL_PARAMETERS = %d\n",
(LONG_INT) OPTIONS->SEQUENTIAL_PARAMETERS);
#endif
#endif
fprintf (ptr_asa_out, "OPTIONS->INITIAL_PARAMETER_TEMPERATURE = %g\n",
OPTIONS->INITIAL_PARAMETER_TEMPERATURE);
fprintf (ptr_asa_out, "OPTIONS->RATIO_TEMPERATURE_SCALES = %d\n",
OPTIONS->RATIO_TEMPERATURE_SCALES);
fprintf (ptr_asa_out, "OPTIONS->USER_INITIAL_PARAMETERS_TEMPS = %d\n\n",
OPTIONS->USER_INITIAL_PARAMETERS_TEMPS);
fprintf (ptr_asa_out, "OPTIONS->TESTING_FREQUENCY_MODULUS = %d\n",
OPTIONS->TESTING_FREQUENCY_MODULUS);
fprintf (ptr_asa_out, "OPTIONS->ACTIVATE_REANNEAL = %d\n",
OPTIONS->ACTIVATE_REANNEAL);
fprintf (ptr_asa_out, "OPTIONS->REANNEAL_RESCALE = %g\n",
OPTIONS->REANNEAL_RESCALE);
#if INT_LONG
fprintf (ptr_asa_out, "OPTIONS->MAXIMUM_REANNEAL_INDEX = %ld\n",
(LONG_INT) OPTIONS->MAXIMUM_REANNEAL_INDEX);
#else
fprintf (ptr_asa_out, "OPTIONS->MAXIMUM_REANNEAL_INDEX = %d\n\n",
OPTIONS->MAXIMUM_REANNEAL_INDEX);
#endif
fprintf (ptr_asa_out, "OPTIONS->DELTA_X = %g\n",
OPTIONS->DELTA_X);
fprintf (ptr_asa_out, "OPTIONS->DELTA_PARAMETERS = %d\n",
OPTIONS->DELTA_PARAMETERS);
fprintf (ptr_asa_out, "OPTIONS->USER_TANGENTS = %d\n",
OPTIONS->USER_TANGENTS);
fprintf (ptr_asa_out, "OPTIONS->CURVATURE_0 = %d\n\n",
OPTIONS->CURVATURE_0);
fprintf (ptr_asa_out, "OPTIONS->QUENCH_PARAMETERS = %d\n",
OPTIONS->QUENCH_PARAMETERS);
fprintf (ptr_asa_out, "OPTIONS->QUENCH_COST = %d\n\n",
OPTIONS->QUENCH_COST);
fprintf (ptr_asa_out, "\n");
}
#if TIME_CALC
/***********************************************************************
* print_time
* This calculates the time and runtime and prints it.
***********************************************************************/
#if HAVE_ANSI
void
print_time (char *message, FILE * ptr_asa_out)
#else
void
print_time (message, ptr_asa_out)
char *message;
FILE *ptr_asa_out;
#endif /* HAVE_ANSI */
{
int who = RUSAGE_SELF; /* Check our own time */
struct rusage usage;
/* get the resource usage information */
#if TIME_STD
syscall (SYS_GETRUSAGE, who, &usage);
#else
getrusage (who, &usage);
#endif
/* print the usage time in reasonable form */
aux_print_time (&usage.ru_utime, message, ptr_asa_out);
}
/***********************************************************************
* aux_print_time
* auxiliary print the time routine
***********************************************************************/
#if HAVE_ANSI
void
aux_print_time (struct timeval *time, char *message,
FILE * ptr_asa_out)
#else
void
aux_print_time (time, message, ptr_asa_out)
struct timeval *time;
char *message;
FILE *ptr_asa_out;
#endif /* HAVE_ANSI */
{
static double sx;
double us, s, m, h;
double ds, dm, dh;
/* calculate the new microseconds, seconds, minutes, hours
and the differences since the last call */
us = (double) ((int) ((double) EPS_DOUBLE + time->tv_usec)) / 1.E6;
s = (double) ((int) ((double) EPS_DOUBLE + time->tv_sec)) + us;
ds = s - sx;
sx = s;
h = (int) ((double) EPS_DOUBLE + s / 3600.);
m = (int) ((double) EPS_DOUBLE + s / 60.) - 60. * h;
s -= (3600. * h + 60. * m);
dh = (int) ((double) EPS_DOUBLE + ds / 3600.);
dm = (int) ((double) EPS_DOUBLE + ds / 60.) - 60. * dh;
ds -= (3600. * dh + 60. * dm);
/* print the statistics */
fprintf (ptr_asa_out,
"%s:time: %gh %gm %gs; incr: %gh %gm %gs\n",
message, h, m, s, dh, dm, ds);
}
#endif /* TIME_CALC */
#endif /* ASA_PRINT */